function out   = doMacroSpanRegressions(allYields,data,setupFilter,numPCAs)

[n1,n2] = size(allYields);
if n2 > n1
    allYields = allYields';
end
[n1,n2] = size(data);
if n2 > n1
    data = data';
end

% Macro variables
T      = size(data,1);
labor  = data(:,1);
wage   = data(:,2);
dc     = data(:,3);
infl   = data(:,4);
r      = data(:,5);
slope  = allYields(:,40)-allYields(:,4);    %long rate minus one-year rate

% Using all bond yields 
%[pcaLoadings, pcaFactors,~,~,pcaR2] = pca(allYields,'NumComponents',6);
% Using only the six bond yields in the estimation
[pcaLoadings, pcaFactors,~,~,pcaR2] = pca(data(:,5:10),'NumComponents',numPCAs);

% Regression I
% M_t = beta0 + beta1'*pca_t + u_t
numBoot = 10000;
window  = 50;
tmp     = nwest(labor,[ones(T,1),pcaFactors(:,1:numPCAs)],4);
if T < 1000
    [R2_se,R2_fractile] = getR2SEboot(labor,[ones(T,1),pcaFactors(:,1:numPCAs)],window,numBoot);
    labor_regI = [tmp.rsqr R2_fractile]*100;
else
    labor_regI = [tmp.rsqr]*100;
end

tmp     = nwest(wage,[ones(T,1),pcaFactors(:,1:numPCAs)],4);
if T < 1000
    [R2_se,R2_fractile] = getR2SEboot(wage,[ones(T,1),pcaFactors(:,1:numPCAs)],window,numBoot);
    wage_regI = [tmp.rsqr R2_fractile]*100;
else
    wage_regI = [tmp.rsqr]*100;
end

tmp    = nwest(dc,[ones(T,1),pcaFactors(:,1:numPCAs)],4);
if T < 1000
    [R2_se,R2_fractile] = getR2SEboot(dc,[ones(T,1),pcaFactors(:,1:numPCAs)],window,numBoot);
    dc_regI    = [tmp.rsqr R2_fractile]*100;
else
    dc_regI    = [tmp.rsqr]*100;
end

tmp        = nwest(infl,[ones(T,1),pcaFactors(:,1:numPCAs)],4);
if T < 1000
    [R2_se,R2_fractile] = getR2SEboot(infl,[ones(T,1),pcaFactors(:,1:numPCAs)],window,numBoot);
    infl_regI  = [tmp.rsqr R2_fractile]*100;
else
    infl_regI  = [tmp.rsqr*100 ];
end

% r_t = beta0 + beta1'*pca_t + beta2'*M_t
%tmp     = nwest(r_t(1:end,1),[ones(T,1),labor(1:end,1),dc(1:end,1),infl(1:end,1)],m); 
% Labor is insignificant in the data!
%tmp      = nwest(r(1:end,1),[ones(T,1),infl(1:end,1),dc(1:end,1)],setupFilter.CSregress_m);
%r_regI   = [tmp.rbar tmp.beta'];


% Regression II
% rx_t+m,k = beta0 + beta1'*pca_t + beta2'*M_t + u_t+m
%where  rx_t+m,k = -(k-m)/4*y_t+m,k-m + k/4*y_t_k - 1/4*y_t_m
m          = setupFilter.CSregress_m;
matSelect  = m:m:max(setupFilter.maturities);
shorRate   = allYields(1:T-m,m);
yieldsUsed = allYields(:,matSelect);
numYields  = size(yieldsUsed,2);
version    = 2;
if version == 1
    % as in Bauer and Rudebusch (2017)
    Y          = zeros(T-m,numYields-1);
    for i=1:numYields-1
        Y(:,i)  = (-m*i/4*yieldsUsed(1+m:T,i) + m*(i+1)/4*yieldsUsed(1:T-m,i+1) - m/4*shorRate);
    end
    X  = [ones(T-m,1),pcaFactors(1:end-m,1:numPCAs),labor(1:T-m,1),wage(1:T-m,1),dc(1:T-m,1),infl(1:T-m,1)];
    %X  = [ones(T-m,1),pcaFactors(1:end-m,1:numPCAs)];
    tmp = nwest(mean(Y,2),X,1);
    xhr_regII          = zeros(1,1+3);
    xhr_regII(1,1)     = tmp.rsqr;
    %xhr_regII(1,2:end) = tmp.beta';
    %tmp.tstat
else
    xhr_regII  = zeros(1,numYields-1);
    for i=1:numYields-1
        % Note that m*i gives k-m in relation to the notation above
        Y  = (-m*i/4*yieldsUsed(1+m:T,i) + m*(i+1)/4*yieldsUsed(1:T-m,i+1) - m/4*shorRate);
        X  = [ones(T-m,1),pcaFactors(1:end-m,1:numPCAs),labor(1:end-m,1),wage(1:end-m,1),dc(1:end-m,1),infl(1:end-m,1)];
        %X  = [ones(T-m,1),pcaFactors(1:end-m,1:numPCAs),dc(1:end-m,1),infl(1:end-m,1)];
        %X  = [ones(T-m,1),pcaFactors(1:end-m,1:numPCAs)];
        tmp = nwest(Y,X,2*4);
        xhr_regII(1,i)     = tmp.rsqr;
        
        %xhr_regII(i,2:end) = tmp.beta(1+numPCAs+1:end,1);
    end
end

% Regression III
% M_t+1 = beta0 + beta1'*pca_t + u_t+1
tmp          = nwest(labor(2:end,1),[ones(T-1,1),pcaFactors(1:end-1,1:numPCAs)],0);
labor_regIII = [tmp.rsqr*100];
tmp          = nwest(wage(2:end,1),[ones(T-1,1),pcaFactors(1:end-1,1:numPCAs)],0);
wage_regIII  = [tmp.rsqr*100];
tmp          = nwest(dc(2:end,1),[ones(T-1,1),pcaFactors(1:end-1,1:numPCAs)],0);
dc_regIII    = [tmp.rsqr*100 ];
tmp          = nwest(infl(2:end,1),[ones(T-1,1),pcaFactors(1:end-1,1:numPCAs)],0);
infl_regIII  = [tmp.rsqr*100];

% M_t+1 = beta0 + beta1'*pca_t + beta2'*M_t + u_t+1
% tmp          = nwest(labor(2:end,1),[ones(T-1,1),pcaFactors(1:end-1,1:numPCAs),labor(1:end-1,1),dc(1:end-1,1),infl(1:end-1,1)],0);
% labor_regIII = [tmp.rsqr tmp.beta(1+numPCAs+1:end,1)'];
% tmp          = nwest(dc(2:end,1),[ones(T-1,1),pcaFactors(1:end-1,1:numPCAs),labor(1:end-1,1),dc(1:end-1,1),infl(1:end-1,1)],0);
% dc_regIII    = [tmp.rsqr tmp.beta(1+numPCAs+1:end,1)'];
% tmp          = nwest(infl(2:end,1),[ones(T-1,1),pcaFactors(1:end-1,1:numPCAs),labor(1:end-1,1),dc(1:end-1,1),infl(1:end-1,1)],0);
% infl_regIII  = [tmp.rsqr tmp.beta(1+numPCAs+1:end,1)'];

%% Correlation of slope with infl and r
tmp         = nwest(slope,[ones(T,1),infl],4);
slopeOnInfl = [tmp.rsqr tmp.beta(2,1) corr(slope,infl)];
tmp         = nwest(slope,[ones(T,1),r],4);
slopeOnR    = [tmp.rsqr tmp.beta(2,1) corr(slope,r)];

%% Ability of slope to forecast inflation and consumption
nn = 4;
%tmp = nwest(dc(1+nn:T),[ones(T-nn,1)  slope(1:T-nn,1)],nn);
tmp = nwest(infl(1+nn:T),[ones(T-nn,1)  slope(1:T-nn,1)],nn);
inflOnSlope = [tmp.rsqr tmp.beta(:,1)' tmp.tstat'];

%dc_yoy    = (dc(1:T-3)+dc(2:T-2)+dc(3:T-1)+dc(4:T))/4;
%infl_yoy  = (infl(1:T-3)+infl(2:T-2)+infl(3:T-1)+infl(4:T))/4;
%slope_ForYoY = slope(4:T);
%Tyoy      = T-3;
%tmp = nwest(dc_yoy(1+nn:Tyoy),[ones(Tyoy-nn,1) slope_ForYoY(1:Tyoy-nn,1)],nn);
%tmp = nwest(infl_yoy(1+nn:Tyoy),[ones(Tyoy-nn,1) slope_ForYoY(1:Tyoy-nn,1)],nn);
%inflOnSlope = [tmp.rsqr tmp.beta(:,1)' tmp.tstat'];

% if T < 1000
%     % Empirical data
%     cDevTrend = setupFilter.cDevTrend/std(setupFilter.cDevTrend);
%     slope_trend = slope;
%     Ttrend = T;
% else
%     % Simulated data
%     if nargin == 5
%         % Use the true value of c from trend
%         cDevTrend = ySim(1,:)'./std(ySim(1,:));
%         slope_trend = slope;
%         Ttrend = T;
%     else
%         % estimate the value of c from trend
%         c = cumsum(dc);
%         resOut = runHamiltonFilter(c);
%         cDevTrend = resOut(13:end);
%         slope_trend = slope(13:end);
%         Ttrend = T-12;
%     end
% end
% tmp = nwest(cDevTrend(1+nn:Ttrend),[ones(Ttrend-nn,1) ...
%             cDevTrend(1:Ttrend-nn,1) slope_trend(1:Ttrend-nn,1) ],nn);
% cOnSlope = [tmp.rsqr tmp.beta(:,1)' tmp.tstat'];
%% Output
out.pcaR2        = pcaR2(1:numPCAs,1)';
out.labor_regI   = labor_regI;
out.wage_regI    = wage_regI;
out.dc_regI      = dc_regI;
out.infl_regI    = infl_regI;
%out.r_regI       = r_regI;
out.xhr_regII    = xhr_regII;
out.labor_regIII = labor_regIII;
out.wage_regIII  = wage_regIII;
out.dc_regIII    = dc_regIII;
out.infl_regIII  = infl_regIII;
out.slopeOnInfl  = slopeOnInfl;
out.slopeOnR     = slopeOnR;
%out.cOnSlope     = cOnSlope;
out.inflOnSlope  = inflOnSlope;
end